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ABSTRACT 


This  thesis  is  concerned  with  a  comparative  study  of  discrete 
time  filters  using  the  theories  of  Wiener-Kolmogorov,  Bode-Shannon, 
and  Kalman,  applied  to  the  filtering  of  a  non-stationary  random 
signal  in  the  presence  of  measurement  noise.   Programs  are  developed 
for  the  simulation  of  these  systems  and  signals  on  a  digital  computer 
Their  filtering  properties  are  compared  for  a  random  input  signal 
with  known  steady  state  characteristics,  starting  at  initial  time 

t  =  t  .   The  results  show  that  when  the  gain  of  the  Kalman  filter 

o  ° 

is  equal  to  its  steady  state  value,  the  Kalman  and  Wiener-Kolmogorov 
filters  perform  identically.   For  large  initial  errors  the  Kalman 
filter,  with  large  initial  gain,  gives  the  best  transient  response; 
for  small  initial  errors  the  Kalman  and  Wiener-Kolmogorov  filters 
are  essentially  equivalent  in  their  transient  responses. 
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CHAPTER  I 

INTRODUCTION 

1.   OBJECTIVE 

This  thesis  examines  and  compares  several  techniques  for  the 
filtering  of  random  signals  in  the  presence  of  measurement  noise. 
In  diagram  form,  the  problem  is  as  shown  in  Figure  1;  x(t)  is  the 
random  signal  to  be  estimated,  v(t)  is  a  random  measurement  noise, 
x(t)  is  the  estimated  or  filtered  signal,  and  e(t)  is  the  error  in 
estimation.   The  criteria  for  the  optimum  filter  is  one  which 
reduces  the  mean  square  error  to  a  minimum.   The  following  techniques 
are  studied  and  compared  with  a  common  example: 

(1)  Discrete  Wiener-Kolmogorov   filter 

(2)  Discrete  Kalman  filter 

(3)  Bode-Shannon  discrete  time  filter 

For  all  of  the  filters  it  is  assumed  that  the  statistical 
properties  of  the  signal  and  the  noise  are  known  by  their  correlation 
functions,  or  their  power  density  spectra,  or  by  equivalent  white 
noise  excited  dynamic  models.   It  is  also  assumed  that  the  signal 
and  the  measurement  noise  are  uncorrelated.   Wiener-Kolmogorov  theory 
assumes  that  the  signal  and  the  noise  are  stationary;  that  is,  their 
statistical  properties  do  not  vary  with  time.   This  implies  that  all 
signals  are  initiated  at  t  =  -».   Kalman  and  Bode-Shannon  filtering, 
on  the  other  hand,  assume  that  signals  start  at  some  initial  time, 
t  ,  and  that  as  time  progresses  the  stationary  probability  charac- 
teristics of  the  signals  are  established.   In  addition,  Kalman  (as 
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well  as  Bode- Shannon)  assumes  that  the  random  signal  process  is 

Markovian,  that  is,  can  be  thought  of  as  generated  by  a  linear 

* 

dynamic  system  excited  by  white  noise.    During  the  transient 

stage,  from  t  =  t„  until  the  statistical  properties  are  established, 
the  Kalman  and  Bode-Shannon  filters  behave  as  time-varying  devices. 
In  the  steady  state,  all  of  the  filters  can  be  expected  to  behave 
identically  since  they  each  establish  a  minimum  mean  square  error 
estimate . 

In  order  to  test  and  compare  the  operation  of  these  filters, 
particularly  for  a  non- stationary  input,  the  Wiener-Kolmogorov  and 
Kalman  filters  have  been  designed  for  discrete  time  simulation  on 
a  digital  computer  for  filtering  the  same  random  input  signal  and 
measurement  noise.   The  output  of  each  filter  and  the  actual  random 
input  signal  have  been  compared  starting  at  time  t  =  t_  until  steady 
state  conditions  are  established. 

The  type  of  filtering  problem  studied  can  be  compared  with  a 
basic  tracking  problem  where  a  target,  with  statiscally  known 
maneuvering  properties,  suddenly  appears.   It  is  desired  to  track 
the  target  and,  in  particular,  to  investigate  how  long  it  takes  the 
output  to  reach  the  steady  state  minimum  mean  squared  error.   The 
results  of  this  investigation  are  presented  in  this  thesis  together 
with  a  summary  of  pertinent  background  theory.   In  addition,  the 
discrete  time  formulations  for  the  Wiener-Kolmogorov  theory  have 
been  developed  and  programmed  for  the  specific  example  treated. 


*R .  E .  Ka Irian ,  New  Methods  and  Results  in  Linear  Prediction 
and  Estimation  Theory,  Technical  Report  No.  61-1,  (Baltimore, 
Maryland:   Research  Institute  for  Advanced  Study,  1961)  p.  3. 
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The  basic  theory  for  discrete  Kalman  and  Bode-Shannon  filtering 
appears  in  the  literature  and  has  been  adapted  to  the  example. 
Programs  necessary  for  generating  the  stochastic  signals  from 
white  noise  have  also  been  developed. 

2.   HISTORICAL  BACKGROUND 

The  many  present  day  theories  for  the  filtering  of  time 
sequences  are  basically  derived  from  the  original  work  of  Wiener 

0   "I  /   1  "7 

and  Kolmogorov  published  in  1941  and  1942.  '   '    They  obtained 
the  direct  solution  of  the  filter  problem  independently  and 
simultaneously  by  the  use  of  the  calculus  of  variations.   This 
approach  led  to  the  filter's  impulse  response  being  characterized 
by  the  solution  of  a  difficult  integral  equation  known  as  the 
Wiener-Hopf  equation. 

One  of  the  next  major  developments  in  the  solution  of  the 

2  17 
filtering  problem  was  accomplished  by  Bode  and  Shannon  :    in  1950 

with  the  simplification  of  the  derivation  of  the  Wiener-Kolmogorov 

filter.   They  used  frequency  domain  analysis  and  conventional 

circuit  theory  concepts  to  interpret  mathematical  operations  in 

physically  intuitive  terms.   Bode  and  Shannon  also  applied  this 

frequency  domain  technique  to  discrete  signals  with  constant 

sampling  intervals. 

In  most  cases,  the  solution  of  the  Wiener-Hopf  integral  equation 

is  a  formidable  task.  With  the  advent  of  the  digital  computer,  an 

alternate  approach  was  sought  to  avoid  this  difficulty.   In  1960 

9  10  11        11 
Kalman  '   '   and  Bucy   attacked  the  filtering  problem  from  the 

state  variable  point  of  view.   They  evolved  a  set  of  differential 


12 


equations,  equivalent  to  the  Wiener-Hopf  integral  equation,  whose 
solution  yields  the  optimum  filter.   These  differential  equations 
may  be  synthesized  in  a  sequential  fashion  and  therefore  digital 
computer  solutions  easily  obtained. 


13 


CHAPTER  II 


THE  WIENER- KOLMOGOROV  FILTER 


There  are  several  methods  for  deriving  the  continuous  Wiener- 
Kolmogorov  filter.  Appendix  A  summarizes  the  solution  which  leads 
to  the  Wiener-Hopf  integral  equation. 


■/ 


cp   (t)  =      H(T)  cp   (t  -  T)dT,   t  £  0  (2.1) 

zx  zz 


H(t)  is  the  desired  filter  impulse  response,  cp   (t)  is  the  auto- 

z  z 

correlation  function  of  the  signal  plus  noise,  and  cp   (t)  is  the 
crosscorrelation  function  of  the  signal  plus  noise,  and  the  signal. 

The  requirements  that  must  be  fulfilled  to  use  Eq.  2.1  are 
(1)  that  all  signals  are  stationary  and  exist  for  all  time;   (2)  that 
the  signal  and  the  noise  are  uncorrelated.   In  order  to  solve  the 
equation,  the  autocorrelation  function  of  the  signal  and  the  auto- 
correlation function  of  the  noise  must  be  known.  These  may  be 
obtained  by  taking  the  Fourier  transform  of  the  respective  power 
spectral  densities,  as  given  by  the  Wiener-Khintchine  equations: 


wx(u)>=  h    (\x^~im^-  k? 


-00 


cp   (T) 
TXX 


George  C.  Newton,  Jr.,  Leonard  A.  Gould,  and  James  F.  Kaiser, 
Analytical  Design  of  Linear  Feedback  Controls  (New  York:   John  Wiley 
and  Sons,  Inc.,  1957)  p.  144. 

John  L.  Stewart,  Fundamentals  of  Signal  Theory  (New  York: 
McGraw-Hill  Book  Company,  Inc.,  1960)  p.  306. 
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cp      (t)   =      I     W    (uo)e   J  T   duo   =  2njf 

XX  X 


w   (oj) 


where  W    (u>)    is   the   power  density   spectrum  of  x(t).      Even  when 

X 

sufficient  knowledge  of  the  signal  is  available,  the  solution  of 
Eq.  2.1  is  generally  an  extremely  difficult  task. 

1.   SOLUTION  BY  SPECTRUM  FACTORIZATION 

The  solution  of  Eq.  2.1  using  the  technique  known  as  spectrum 

factorization  has  been  developed  by  several  authors.  '   '    The 

autocorrelation  function  of  the  signal  plus  noise,  cp   (s),  is 

zz 

factored  so  that 


CPZZ(S)  =CP  ZZ(S)CP  zz(8) 


(2.1-1) 


where  cp   (s)  contains  those  factors  of  cp   (s)  which  have  poles  and 
T  zz  zz 

zeros  in  the  left  half  plane.  The  term 


CPZX(S) 


(2.1-2a) 


+ 


is  defined  by  those  terms  of  the  partial  fraction  expansion  of: 


cp   (s) 

Tzx 


*  zz(s) 


which  have  poles  in  the  left  half  of  the  s-plane.  That  is 


CP 


zx 


9 


zz 


9   (s) 

zx 

zz 


zx 


(2.1-2b) 
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The  resulting  weighting  or  transfer  function  of  the  filter  is  given 

by: 

cp   (s) 

zx 


H(s)  = 


zz 


+ 


9   (s) 
T  zx 


(2.1-3) 


The  use  of  this  equation  may  be  demonstrated  by  the  example 
shown  in  Figure  1.   Consider  that  x(t)  has  an  autocorrelation 
function: 


9xx(t)  =  2  e 


(2.1-4) 


The  Laplace  transform  of  the  autocorrelation  function  is  then 


cp   (s)  = 
Txx 


1  -  s 


(2.1-5) 


which  yields  the  power  density  spectrum, 


wxw  ■  *r 


1  +  10 


(2.1-6) 


This  signal  can  be  shown  to  correspond  to  the  random 

walk  signal  shown  in  Figure  2,  where  x(t)  is  a  rectangular  wave  with 

values  +  —  or  -  —  and  with  zero  crossings  located  at  event  points 

which  are  Poisson  distributed  with  an  average  frequency  of  —  crossings, 

per  second.   This  time  domain  signal  is  not  unique  since  many  signals 

may  have  the  same  power  density  spectrum. 

Now  consider  a  white  noise  source  where, 

(Pw(t)  =  6(t)  (2.1-7) 

where  6(t)  is  the  Dirac  delta  function,  which  has  the  following  Laplace 
transform: 

<Pvv(s)  -  1  (2.1-8) 

16 
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The  autocorrelation   function  of   the   signal   plus   noise   is 


cp      (s)    =  cp      (s)   +  cp      (s)   +cp      (s)    +  cp      (s)  (2.1-9) 

ZZ  XX  Tw  Tvx  Txv 


Since  the  signal  and  noise  are  assumed  to  be  uncorrelated  the  last 

two  terms  are  zero  and  cp   (s)  reduces  to: 

Tzz 

cp   (s)  =  cp   (s)  +  cp   (s) 
Tzz      Txx      Tvv 


1  2  -  s2 

+  1  = ~-  (2.1-10) 


2  i     2 

1  -  s  1  -  s 


The  crosscorrelation  function  of  the  signal  plus  noise,  and  the  signal 
is  : 

00 

Cp   (T)    =      [   (V(t)  +  X(t))   X(t+T)dt 
ZX  J 

^00 

00  00 

v(t)  x(t  "-T)dt         x(t)  X(t  +  T)dt 


I  v(t)  X(t  +  T)dt   +    I 


x(t)  x(t  +  T)dt   =  cp   (t)  (2.1-11) 

XX 


Since  the  noise  and  the  signal  are  uncorrelated, 

1  -  s 
Factoring  Eq.  2.1-10  results  in: 

cp   (s)  -  ^d-  ■   (V7+S)   (V7°  S)  (2-1"13) 

ZZ       1  -  s2        (1  +  s)     (1  -  s) 


Therefore,  from  Eq.  2.1-1, 
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and, 


zz 


V7 


+  s 


1  +   s 


*  zz(s) 


VT- 


1  -    s 


Now, 


CP«x(s) 
*"«z(8) 


1  1   -    s 

i  -  s2    (\TT- 


s)  (1  +  s)      (  V2   -    s) 


+ 


2.414         WpT *    s  1  +  s 


Therefore,  fromEq.  2.1-2a, 


(2.1-14) 


CPzx(s) 
cp"zz(s) 


1 


2.414     1  +  s 


(2.1-15) 


Substituting  into  Eq.  2.1-3, 


H(s)    = 


1 


2.414 


1       I     I    1  +  s 
1  +  s        \j2  +  s 


2.414       k/2~+  s 


(2.1-16) 


The  resulting  filter  weighting  function,  or  impulse  response,  is 


1     ,-^,  tM 


H 


(t>- 


2.414 


0 


,   t  <  0 
This  filter  is  physically  realizable,  since  it  is  causal. 


(2.1-17) 
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THE  BODE- SHANNON  SOLUTION 


The  frequency  domain  technique  for  the  Wiener-Kolmogorov  filter 
as  developed  by  Bode  and  Shannon,   is  to  convert  the  signal  power 
spectral  density  to  that  of  white  noise,  and  then  have  a  filter 

operate  on  this  white  noise .   The  development  of  this  technique 

2,14,17    ,  .   ,.       .  .   . 
appears  in  several  references         and  is  discussed  in  Appendix 

B.   The  following  results  are  obtained. 

Given  W  (w),  the  power  spectral  density  of  the  signal  x(t) ; 

X 

and  W  (uo)  ,  the  power  spectral  density  of  the  noise  v(t)  ,   consider 
the  transfer  function  of  the  following  minimum  phase  (zeros  and 
poles  in  the  left  half  of  the  s-plane)  filter: 

B(u>)  =   ^/W  (u>)  +  W  (u3)  (2.2-1) 

From  a  white  noise  input  this  filter  will  produce  a  power  spectral 
density  equal  to  that  of  the  signal  plus  noise.   As  shown  in  Appendix 
B,  the  following  transfer  function  can  then  be  obtained  for  the  non- 
causal  optimal  filter. 

W  (w)  B(u)) 

h  (u>)  =   r!!  />,,„/  m  (2.2-2) 

uo       [W  (a))  +  W  (cjo)] 

X  V 


The  preceding  equations  are  represented  in  Figure  3.   Since  h 0*0  is 
non-causal,  let 


Ralph  J.  Schwarz   and  Bernard  Fried  land,  Linear  Systems 
(New  York:   McGraw-Hill  Book  Company,  1965),  p.  334-45. 
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%>M  = 


,  t  <  0 


h  (t)   ,  t  *  0 

uu      ' 


(2.2-3) 


where  g  (t)  represents  the  portion  of  h  (t)  which  exists  for  t  ^  0. 
This  is  equivalent  to  the  spectral  factorization  of  Eq.  2.1-1  and 
Eq.  2.1- 2a.   The  desired  optimum  causal  filter  has  the  transfer 
function 


H(u>)  = 


G  (uu) 
B  (u>) 


(2.2-4) 


Now  consider  the  example  problem  again.   The  signal,  x(t),  has 
a  power  spectral  density 


W  (ou)  = 


_1_ 

2rr 


1  +  uu' 


(2.2-5) 


and  the  uncorrelated  measurement  noise  is  white,  with  a  power  spectral 
density 


W  (Oi)  =  ~ 

v      2n 


(2.2-6) 


The  power  spectral  density  of  the  signal  plus  noise  is  therefore 


W  (uu)  =  W  (uu)  +  W  (uu)  =  — 

2TT 


This  may  be  factored 


+  1 


I 


1  +  uu' 


l_ 

2u 


2  +  uu 
1  +  m 


(2.2-7) 


zk  ;    2tt    .(1  +  juo)     (1  -  joj) 


(2.2-8) 


Thus,  the  minimum  phase  transfer  function  required  to  produce  this 
spectrum  from  white  noise,  using  Eq.  2.2-1,  is: 
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vr 


+  jou 


B(u))  -  m  ^rt  (2-2-9) 

The  inverse  of  this  transfer  function  will  produce  white  noise  from 
the  signal  plus  noise. 

Now,  substituting  in  Eq .  2.2-2, 

h    (u))      =       2tt  \i  +u)2/\    i  +   yu/      V^ 

ur   '  2 

1  2  +  u) 


2 

2tt  1  +  uo 


sfl   + 


JHL 


(2  +  uo2)    (1  +  ju>)         n/Itt 


CfT-     p)     (1    +    jU>)     V2TT 


1  +  1        \  (2.2-10) 

2.414  V2rr     I    y/T  -    ju>  1  +  ju) 


By   taking   the   inverse   Fourier   transform  of  Eq.    2.2-10,    the   impulse 
response,    or  weighting   function,    is   then 

1  J  ?      t 

I  P=         — - —       e_t  ,    t  *  0  (2.2-11) 

|\/2tt  2.414 

( 

which  is  non-causal.   The  best  causal  impulse  response,  from  Eq .  2.2-3, 
is 
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,    t  <   0 


V2tt        2.414     e 


,    t  *  0 


The   desired   impulse  response   is    therefore, 


H(uo)    = 


G   (u>) 
B(«u) 


!       1         'i 


/Ztt     \ 2.414  I     I  1  +  p 


I 


tff\  ^ 


(2.2-12) 


and 


2.414       /F+  jo> 


H(t)    =, 


V2 


2.414 


,    t  <  0 


,   t  *  0 


(2.2-13) 


This    filter    impulse   response    is    identical    to  Eq .    2.1-17,    as    expected. 

3.      DISCRETE   TIME    SOLUTION  OF  THE 
WIENER- KOLMOGOROV  FILTER 

The    initial   problem   to  be    solved   for   digital   computer   simulation 
of   the  Wiener- Kolmogorov   filter    is   the   realization   of   the  discrete 
signal,    x(kT),    based   on   the   required   power    spectral   density 


W    (oj)    =     — 

2TT 


b 


(2.3-1) 


1+0) 

Since  the  power  density  spectrum  of  white  noise  is  a  constant,  this 
power  spectral  density  may  be  obtained  by  driving  a  linear  system, 
G,(juO,  by  white  noise,  w(t)  ,  so  that 
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W  ((JO)  =  GAp)    GA-P)   W  (uo) 
X        1       1       w 


2tt 


l+O) 


J 


(2.3-2) 


where 


GAP)    = 


1  +  p 


It  is  extremely  difficult,  if  not  impossible,  to  obtain  a 
discrete  white  noise  source  which  is  flat  over  all  frequency.   In 
practice,  what  needs  to  be  done  is  to  make  the  spectrum  of  the  noise 
source  flat  over  the  bandwidth  of  G,  ( p)  .   The  3  decibel  cutoff 
frequency  for  G1 ( p)    is  one  radian  per  second.   The  effective  noise 
bandpass  is  given  by 


B  =  —  (1)  =  —  radian/sec.  =  1/4  Hz. 


(2.3-3) 


TT  * 

where  noise  bandwidth  is  —  times  the  3  decibel  cutoff  frequency.  . 
The  sampling  rate  for  the  discrete  noise  should  be  much  greater  than 
1/4  Hz.,  and  100  Hz.  was  chosen  as  convenient.   The  sampling  period, 

T,  therefore  is  0.01  seconds.   The  spectrum  of  the  noise  should  then 
be  flat  to  50  Hz.  as  shown  in  Figure  4. 

The  autocorrelation  function  of  this  band  limited  white  noise 
is  then  given  by 


,  .     100  sin  (100  ttt) 
^(T)  =  100^ 


WW 


(2.3-4) 

as  shown  in  Figure  5.   The  variance  of  w(t),  the  noise  sequence,  is 
given  by 

<P   (0)  =  o-2  =  100  =  |  (2.3-5) 

ww       w         T 


Mischa  Schwartz,  Information  Transmission,  Modulation,  and 
Noise  (New  York:  McGraw-Hill  Book  Company,  Inc.,  1959),  pp.  206-207 
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AUTOCORRELATION    OF    WHITE    NOISE   SOURCE 
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The  white  noise  source  is  then  simulated  by  a  series  of  pulses  of 
repetition  period  .01  seconds,  whose  amplitudes  have  a  Gaussian 
distribution  with  zero  mean  and  variance  100,  followed  by  a  zero 
order  hold  as  shown  in  Figure  6.   The  recursive  equation  necessary 
for  the  computer  simulation  of  the  desired  signal  source  is  derived 
by  taking  the  Z  transform  of  the  input-output  transfer  function  of 
Figure  6. 


x(z) 
w(z) 


=  Z 


s  +  1 


sT 


(2.3-6) 


Hence 


and 


-T  -T 

x(z)  (z  -  e   )  =  w(z)  (1  -  e   ) 


-T  -T 

zx(z)  =  x(z)e   +  w(z)  (1  -  e   ) 


Therefore, 


x(k  +  1)  =  e"Tx(k)  +  (1  -  e"T)  w(k);  T  =  .01  sec.     (2.3-7) 


To  check  that  the  signal  does,  in  fact,  have  the  desired  power 
spectral  density,  a  subroutine  entitled  HARM  was  used  to  obtain  the 
Fast  Fourier  Transform.   The  Fourier  coefficients  were  then  squared 
and  plotted  versus  radian  frequency,  and  the  graph,  Figure  7,  obtained 
Included  on  this  graph  is  the  desired  power  spectral  density  curve. 


International  Business  Machines  Corporation,  System  /360 
Scientific  Subroutine  Package,  (360A-CM-03X)  Version  II  Programmer's 
Manua 1 ,  (White  Plains,  New  York:  International  Business  Machines 
Corporation,  1967) 
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The  program  was  run  for  2048  samples  and  it  can  be  seen  that  the 
power  density  spectrum  approximates  the  desired  spectrum.   The 
computer  program  for  generating  the  signal  and  checking  its  power 
density  spectrum  is  included  as  Appendix  C. 

Now  that  the  proper  signal  has  been  obtained,  a  recursive 
equation  for  the  filter  is  needed.   The  filtering  system  is  drawn 
in  Figure  8.   Using  the  Z  transform,  the  recursive  filter  equation 
is  developed  as  follows: 

-sTl 


s  +  V2" 


(2.3-8) 


Hence 


x(z)   z  -  e 


2  + 


7* 


y  2  T  |  z(z)    (2.3-9) 


and 


z4(z)  =  x(z)  e"  ^  T  +  — ±7=-      (  1  -  e'  ^  T  | 


2  + 


V2 


Therefore,  the  recursive  equation  is: 
X(k  +  1)  =  e"^2  T  2(k)  +  -±+ 


TTJi    V-  e_VrTi  z<k>       <2-3"u) 


(2.3-10) 


The   error   is 


E(k)    =  X(k)    -   X(k) 


(2.3-12) 
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A 

The  assumption  is  made  that  X(0)  is  zero.   The  sampling  interval 
is  0.01  seconds.   The  measurement  noise,  v(t) ,  is  white,  with  a  power 
spectral  density,  W  Qx>)    =  — ,   and  is  developed  in  an  identical  fashion 
to  w(t),  the  driving  noise  for  the  signal  generator.   The  two  noise 
sources  are  uncorrelated.   The  mean  square  error  and  the  variance 
of  the  error  are  then  calculated,  on  a  continuing  basis,  using  the 
equations : 

1  n 

E(n)  =^2    E(k)  (2.3-13) 

k=l 


E(n)2  =-   S    (E(k))2  (2.3-14) 

n   k=l 


a     =  E2(k)  -  E(k)  (2.3-15) 

e 

The  computer  run  was  for  30  seconds,  i.e.,  3000  samples,  so  that 

steady  state  is  reached.   The  results  are  discussed  in  Chapter  V. 

The  computer  program  for  the  Wiener-Kolmogorov  filter  appears  as 

Appendix  D. 
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CHAPTER  III 

BODE- SHANNON  DISCRETE  TIME  FILTERING 

The  Bode  and  Shannon  frequency  domain  approach  may  be  applied 
in  the  time  domain  directly,  as  developed  by  Schwarz  and  Friedland. 
The  formulation  is  now  in  terms  of  summations,  discrete  signals,  and 
discrete  power  spectral  densities.   This  technique  results  in  the 
following  equation  which  is  similar  to  the  Wiener-Hopf  integral 

equation. 

n 
cp   (k,n)  =  2    h(n,j)  cp   (j,k)  (3.1a) 

ZX  *       r\  Z  2 

J=0 


k  =  0,  1,  2,  ... .  n 


where 


cp   (k,n)  =  E   z(k)  x(n) 
zx 


(3.1b) 


is  the  discrete  crosscorrelation  function  between  z(k)  and  x(n), 
and 


cpzz(j,k)  -E 


z(j)  z(k) 


(3.1c) 


is  the  discrete  autocorrelation  function  of  z(k).   The  unit  impulse 
response  of  the  filter  is  h(n,j).   The  notation  h(n,j)  denotes  the 
response  at  time  n  due  to  an  impulse  applied  at  time  j  where  j  ^  n. 
Since  the  impulse  response  of  the  filter  must  be  causal,  this  requires 
that  h(n,j)  equal  zero  for  n  less  than  j.   Equation  3.1  holds  for  both 
stationary  and  non- stationary  signals.   When  the  signals  are  stationary 
then 

cp   (k,n)  =  cp   (n  -  k)  (3.2a) 

Z  X  Za 

<PZZUA)  =cPzz<k  -  j)  <3-2b> 
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h(n,j)  =  h(n  -  j) 


(3.2c) 


and  Eq.  3.1  reduces  to 


cp   (1)  =   2   h(m)  cp   (m  -  l) 

TZX  -  zz 

m=0 


(3. 2d) 


where 


i  =  n  -  k    and  m  =  n  -  j 


(3.2e) 


The  power  spectral  density  of  a  discrete  signal,  z(kT),  is  defined 


as 


W  (z)  =  T  Z 
z 


<f>   (nT) 
zz 


=  T   S   cp   (nT)  z 
n=-»  Tzz 


-n 


(3.3) 


This  is  the  discrete  form  of  one  of  the  Wiener-Khintchine  equations. 
The  prime  denotes  the  power  spectral  density  of  a  discrete  signal. 
It  should  be  noted  that  the  z  variable  within  the  parenthesis  denotes 
the  Z-transform  shifting  variable.   Proceeding  as  in  Appendix  B,  a 

transfer  function,  B(z),  is  developed  which  will  convert  a  signal  with 

i 

a  power  spectral  density  of  W  (z)  to  a  white  noise  process,  u(kT). 

z 


This  requires  that 

B(z)  B(z_1)  w'(z)  =  1 
z 

Since  B(z)  is  causal,  its  Z  transform  must  be  given  by 


(3.4) 


B(z)  =   S   b(n)  z 
n=0 


■n 


(3.5) 


*  17 

Schwarz  and  Friedland   use  this  definition, 

of  1/T  instead  of  T. 


12 
Kuo  "  uses  a  factor 
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If  W  (z)  can  be  factored  into  a  product, 


W  (z)  =  W  (z)  W  (z"1) 


(3.6) 


then 


B(z)  =   - 


W  (a) 
z 


(3.7) 


For  a  white  noise  input  cp   (j,k)  equals  6(j,k),  and  Eq.  3.1  reduces  to 


cp   (k,n)  =  h(n,j)     for  j  i  n 


cp   (k,n)  -  0 
Tux 


for  j  >  n 


(3.8) 


The  impulse  response  of  an  optimal  causal  filter  for  a  white  noise 
input,  u(kT),  and  a  desired  output,  X(kT),  is  then  given  by 


g(n,j)  =  ^(k.n) 


(3.9) 


For  stationary  processes  Eq.  3.9  reduces  to 


g(m)  =  cp   (m)    for  m  =  n  -  k  >  0 

ux 


(3.10) 


Taking  the  Z  transform  of  Eq.  3.10,  with  the  aid  of  Eq.  3.3,  yields 


g(z) 


W   (z) 

ux 


(3.11) 


where  W   is  the  cross  power  spectral  density  of  the  white  noise,  u(kT), 
and  the  signal.   The  subscript  _c  denotes  the  causal  part  of  the  transform 
(terms  involving  z   ) .   However 


W   (z)  =  W   (z)  B(z"1) 
ux       zx 


(3.12) 


so  that 
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g(z)  ■    ' 


Wzx(z)  B(Z"1) 


(3.13) 


The  optimum  filter  is  thus 


H(z)  =  G(z)  B(z) 


w'  (z)  BCz"1) 
zx 


B(z) 


(3.14) 


The  procedure  of  taking  "the  causal  part  of"  is  analogous  to  the 
spectrum  factorization  done  in  Chapter  II. 

Using  the  above  technique,  the  previous  example  is  now  solved 
for  discrete  time  signals.   The  signal  and  noise  correlation  functions 
are 

*xx(kT)  =  2  e"'k'T  (3'15) 


cp  (kT)  =  6(k) 
w 


(3.16) 


The  sampled  power  spectral  density  of  the  signal  may  be  obtained  by 
means  of  Z  transforms  and  is  as  follows: 


W  (x)  =  T  Z 
z 


=  T 


1  -  kT 
26 


1  -  e' 


(1  -  ez  )  (1  -  ez) 


(3.17) 


where  e  equals  e   .  For  the  noise,  the  sampled  power  spectral  density  is 

w'(z)  =  T  (3.18) 
v 

Therefore,  the  power  spectral  density  of  the  signal  plus  noise  input 


to  the  filter  is 


W  (z)  =  w'(z)  +  w'(z) 
z       x       v 


=  T 


-1 
2  -  sz   -  ez 

(1  -  ez"1)  (1  -  ez) 


(3.19) 
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This  may  be  written  in  the  following  manner  for  ease  in  factoring: 


w,(z)=l£   (i  -  bz"1)  (1  -  bz) 
■z>   }  b   (1  -  ez  I)  (1  -  ez) 


(3.20) 


where 


b  = 


1-VT-e2 


(3.21) 


The  transfer  function  needed  to  produce  white  noise  from  the  signal 
and  noise  input  to  the  filter  is  then, 

1 


B(z)  =  Vb/eT 


1  -  ez 


1  -  bz 


-1 


W  (z) 


(3.22) 


Since   the   signal  and   the  noise  are  uncorrelated, 


and 


Cpzx(k)    =  <Pxx<k) 


W      (z)    =  W    (z) 

zx  X 


(3.23) 


(3.24) 


Thus,    from  Eq.    3.13,^ 
g(z)    - 


w'(z)   B(z_1) 

X 


1    -    6' 


(1   -   ez"    )    (1  -   ez) 


\/b/eT 


(1  -  ez) 
(1  -  bz) 


(3.25) 


Jc 


Simplifying  and   taking  the   causal   part  yields 

2 


1  -   e' 


g(z)    =    Vb/eT     ji 


1 


be  1  -   ez 


-1 


(3.26) 


The  discrete  time  transfer  function  for  the  optimum  filter  is  then 
H(z)  =  g(z)  B(z) 


1  -  e 
Te  I  1  -  be 


(3.27) 


1  -  bz 
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The  optimum  filter  diagram  is  shown  in  Figure  9.   Since  b,  e, 
and  T  are  constants,  H(z)  may  be  rewritten  as: 

H(z)  =  C (3.28) 

1  -  bz 

where 

The  recursive  equation  for  the  filter  is  developed  as  follows: 

H(Z)  =  ksl  =       c 

()    Z<2>     l-bz-1 
x(z)  (1  -  bz"1)  =  Cz(z) 

x(z)  =  bz"1  x(z)  +  Cz(z) 

x(k  +  1)  =  bx*(k)  +  CZ(k  +  1)  (3.30) 


The  error  is 


E(k)  =  X(k)  -  X(k)  (3.31) 


For  the  specific  example  the  coefficients  for  the  filter,  Eq.  3.30, 
become   b  =  .869  and  c  =  12.45.   These  differ  from  the  coefficients  for 
the  discrete  Wiener -Kolmogorov  filter,  Eq.  2.3-11.   One  might  expect 
that  these  results  would  be  identical,  but  although  Schwarz  and  Friedland 
consider  both  developments  the  discrepancy  is  not  explained  and  remains 
a  subject  for  further  investigation. 
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CHAPTER  IV 

THE  KALMAN  FILTER 

1.   DISCRETE  TIME  SOLUTION 

The  objective  in  this  chapter  is  to  examine  the  filtering 
problem  from  the  state  point  of  view  as  developed  by  Kalman.  '   ' 
In  the  previous  cases,  the  signal  was  characterized  by  either  the 
autocorrelation  function  or  the  power  spectral  density.   Here  the 
signal  is  characterized  by  the  output  of  a  linear  dynamic  system 
driven  by  white  noise  similar  to  Bode-Shannon.   This  presents  no 
problem,  for  if  the  power  spectral  density  is  rational,  it  can 
always  be  represented  as  the  output  of  a  linear  system  driven  by 
white  noise.   As  in  Chapter  II,  the  linear  system  has  the  transfer 
function, 

Since  this  is  a  first  order  system,  the  observation  matrix  is  unity 
and  scalar  notation  can  be  used.   The  corresponding  differential 
equation  is 

^&-     =  -x(t)  +w(t)  (4.1-2) 

Then,  from  state  variable  theory, 

*  =  e"T    and    T  =  1  -  e"T 

and  the  differential  equation  may  be  expressed  in  equivalent  recursive 
form: 

X(k  +  1)  =  e"T  X(k)  +  (1  -  e"T)W(k)  (4.1-3) 
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The  signal  model  is  outlined  in  diagram  form  in  Figure  10;  w(t) 
and  v(t)  are  independent  gaussian  random  processes  with  zero  means 
and  covariances  given  by: 

cov  [w(t),  w(t)]  =  Q6  (t  -  t) 
cov  [v(t),  v(t)]  =  R6  (t  -  t) 
cov  [w(t),  v(t)]  =  0 

Q  and  R  are  the  variances  of  the  signal  white  noise  source  and  the 
measurement  noise,  which  are  set  equal  to  unity  in  order  to  duplicate 
the  model  for  the  Wiener  filter.   The  noise  sources  are  simulated  in 
the  exact  manner  as  described  in  Chapter  II,  Section  4. 

The  Kalman  filter  '  is  shown  in  Figure  11.   G(k)  is  an  adjustable 
gain  which  minimizes  the  mean  square  error, 

E2(n)  =  -   S    (E(k))2  (4.1-4) 

n  k-1 


where 


E(k)  -  X(k)  -  X(k) 


G(k)  is  determined  by  the  following  equations: 

r(k)   =  P(Vk  -  1) (4  1-5) 

&(k)    P(k/k  -  1)  +R  ^-l   ^ 

P(k/k)  =  (1  -  G(k))  P(k/k  -  1)  (4.1-6) 

P(k  +  1/k)  =  $2  P(k/k)  +  QD  (4.1-7) 


where  P(k/k  -  1)  denotes  P  at  time  k  given  the  value  at  time  k  -  1. 
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P(k/k-l)  =E  I      [z(k)-Z(k/k-l)J  2  I 


1  k~l 
n=0 


Z(n+1)  -  Z(n+l/n) 


P(k/k)  =  E 


Z(k)  -  Z(k/k) 


1   k 

*  -1 

n=l 


Z(n)  -  Z(n/n) 


also, 


qd  =  r  Q 


It  is  necessary  to  define  P(l/0)  in  order  to  start  the  recursive 
process. 

It  is  important  to  note  that  once  P(l/0)  is  specified  the  gain 
equations,  Eqs .  4.1-5,  4.1-6,  and  4.1-7,  are  recursive  and  do  not 
depend  on  the  observations  Z(kT).   They  depend  upon  the  variance 
of  the  signal  driving  noise,  the  variance  of  the  measurement  noise, 
and  the  characteristics  of  the  linear  system  used  for  signal  generation, 

From  Figure  lj,  the  recursive  equation  for  the  Kalman  filter  may 
be  expressed  as 


X(k)  -  $X(k  -  1)  +  G(k) 


Z(k)  -  $X(k  -  1) 


(4.1-8) 


The  error  is 


E(k)  =  X(k)  -  X(k) 


(4.1-9) 


As  in  the  Wiener  simulation,  the  sampling  interval  is  0.01 
seconds,  X(0)  is  the  initial  estimate  at  time  t  =  0  and  the  length 
of  the  computer  run  is  3,000  samples.   The  computer  program  for 
the  Kalman  filter  simulation  is  included  as  Appendix  E. 
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2.   EQUIVALENCE  OF  KALMAN  AND  WIENER- KOLMOG OR OV  FILTERS 

IN  THE  STEADY  STATE 

It  has  been  shown  by  C.  E.  Hutchinson  that  the  continuous 
Kalman  and  Wiener  filters  are  equivalent  in  the  scalar  case  after 
the  steady  state  condition  is  reached,  i.e.,  G(k)  is  a  constant. 

This  equivalence  is  developed  as  follows. 

a2 
Given  a  message  spectrum,  W  (uo)  = 


bV  +  i 


2 
a  noise  spectrum,  W  (w)   ■  c 

v 


and  a  crosscorrelation  function,  9   (t)  =  0 

xv 


The  optimal  causal  filter,  in  the  Wiener  sense,  is  then 


where 


H(s)  =  (k2  -  kx)    s  +  \  (4.2-1) 


kt  =  1/b 


k2  =  1/bc   *  a2  +  c 


2 


Now  consider  the  Kalman  filter.  The  signal  model  is,  in  general, 

defined  by: 

x(t)  =  Fx(t)  +  w(t)  (4.2-2) 

z(t)  =  Mx(t)  +  v(t)  (4.2-3) 

The  best  estimate  of  x,  designated  x\  is  obtained  from  the  continuous 

filter  equations  which  correspond  to  the  discrete  equations  4.1-5, 

4.1-6,  4.1-7,  and  4.1-8. 
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where 


and 


x(t)  =  Fx(t)  +  g(t)  [z(t)  -  Mx(t)] 


g(t)  =  P(t)  mV1 


P(t)  =  FP(t)  +  PfT(t)  -  P(t)M  R_1MP(t)  +  Q 

cov  [w(t),  w(t)]  =  Q6(t  -  t) 
cov  [v(t),  v(t)]  =  R6(t  -  t) 
cov  [w(t),  v(t)]  =  0 

If  a  one  dimensional  problem  is  considered  with 


(4.2-4) 


f 

=  -L/b 

Q 

2/K2 

=  a  /b 

R 

2 
=  c 

M 

=  1 

the 

Kalman  solution 

becomes 

• 

/,   'l/.s 

where 


^(t)  =  -1/b  x(t)  +  g(t)  [z(t)  -  x(t)] 


g(t)  =  P(t)/c^ 

P(t)  =  -  (2/b)P(t)  -  (l/c2)P2(t)  +  a2b2 


(4.2-5) 


In  the  steady  state,  when  P(t)  is  identically  zero,  the  equation  for  the 
Kalman  filter  reduces  to 

x(t)  =  -1/b  x(t)   + 


-1/b  +  1/bc  V/a2  +  c2 


z(t)  -  x(t) 


J 


(4.2-6) 
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A 

The  transfer  function  relating  the  estimated  signal,  X(s),  to 
the  signal  plus  the  noise,  Z(s),  can  be  expressed  as 

which  is  identical  to  Eq.  4.2-1.   This  result  can  be  used  as  a  check 
on  computer  programs  for  the  Kalman  and  Wiener-Kolmogorov  filters. 
When  G(k)  has  reached  a  constant  value,  the  mean  square  errors  should 
be  identical. 
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CHAPTER  V 
RESULTS  AND  COMPARISONS 

The  continuous  Kalman  and  Wiener-Kolmogorov  filters  have  been 
shown,  in  Chapter  IV,  to  be  equivalent  in  the  steady  state  condition, 
that  is,  when  the  Kalman  gain  is  constant.   This  result  is  verified 
experimentally,  as  shown  by  the  computer  curves  of  Figures  16,  20,  24, 
and  28.   The  steady  state  gain  of  the  Kalman  filter  reaches  a  value 
of  .0041335  as  expected.  As  soon  as  the  gain  reaches  its  steady 
value,  the  curves  of  mean  square  error  versus  time  (Figures  18,  19, 
22,  23,  26,  27,  30,  and  32)  are  identical  to  each  other  and  to  the 
Wiener-Kolmogorov  filter  curves  (Figures  14  and  15) . 

In  order  to  start  the  recursive  process,  the  Kalman  filter 
requires  a  priori  knowledge  of  the  initial  error  covariance,  P..  ,„. 
Four  different  values  of  P  ._  (0.0,  .414065,  1.0,  and  10.0)  were 
chosen  to  compare  the  effectiveness  of  the  Kalman  filter  for  different 
values  of  this  term.  The  value  of  P  ,_  =  .415065  was  chosen  because 
it  gives  the  steady  state  gain  (.0041335)  at  the  start  of  the  computer 
run.  With  this  P  /f)  the  performance  of  the  Kalman  and  Wiener-Kolmogorov 
filters  are  identical  (Figures  13-15  and  21-23). 

In  simulating  the  non-stationary  signal,  two  different  initial 
values  were  chosen;  one  to  correspond  to  a  small  initial  error  in 
estimation  and  the  other  to  correspond  to  a  large  error  in  the  initial 
estimation.   In  a  tracking  problem,  the  first  case  would  represent 
the  situation  where  a  target  is  picked  up  at  long  range.   The  second 
case  represents  the  situation  where  a  target  suddenly  appears  at 
relatively  short  range.   For  small  initial  errors,  the  Wiener-Kolmogorov 
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and  the  Kalman  filters  performed  essentially  the  same  (Figures  13, 
14,  17,  18,  21,  22,  25,  26,  29,  and  30).   For  large  initial  errors 
the  Kalman  filter,  with  large  initial  covariance  of  error,  P 1  /n, 
gives  the  best  transient  response  (Fig.  31).   For  large  initial 
error,  if  Pn  ,_  is  equal  to  zero,  the  transient  response  of  the 
Kalman  filter  takes  somewhat  longer  time  to  settle  than  the  Wiener- 
Kolmogorov  filter  as  shown  in  Fig.  19. 
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CHAPTER  VI 
CONCLUSIONS 

The  responses  of  the  discrete  Kalman  filter  for  various  initial 
error  covariance  and  the  discrete  Wiener-Kolmogorov  filter  have  been 
compared  for  a  non-stationary  random  input  with  defined  statistical 
properties.   It  was  necessary  to  develop  a  computer  program  for  the 
simulation  of  the  desired  spectrum  from  a  white  noise  source.   The 
discrete  time  Wiener-Kolmogorov  filter  was  also  programmed  for  the 
digital  computer. 

It  has  been  shown  that  when  the  gain  of  the  Kalman  filter  is 
equal  to  its  steady  state  value  the  performance  of  the  Wiener- 
Kolmogorov  and  Kalman  filters  are  identical.   For  large  initial  error 
the  Kalman  filter  (with  a  large  initial  gain)  shows  the  best  transient 
response.  When  the  initial  error  is  small  the  Wiener-Kolmogorov  and 
the  Kalman  filters  respond  in  as  essentially  equivalent  fashion. 
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APPENDIX  A 

DEVELOPMENT  OF  THE  WIENER-HOPF  INTEGRAL  EQUATION 

The  Wiener-Hopf  equation  gives  the  relationship  between  the 
filter  weighting  function,  the  autocorrelation  function  of  the 
input  to  the  filter,  and  the  crosscorelation  function  of  the  filter 
input  with  the  desired  filter  output.   The  problem  is  described  in 
Figure  1;  x(t)  is   some  signal,  v(t)  is  white  noise,  x(t)  is  the 
estimate  of  the  signal,  and  e(t)  is  the  error  in  estimation;  reduce 
e(t)  to  a  minimum  in  the  mean  squared  sense  by  a  correct  choice  of 
h(t),  the  impulse  response  of  the  filter.   The  error  by  definition  is 


e(t)  =  x(t)  -  x(t) 


(A.l) 


Squaring  both  sides  yields, 


e2(t)  =  x2(t)  -  2x\t)  x(t)  +  x2(t) 


(A. 2) 


The  output,  x(t),  and  the  impulse  response  of  the  filter,  H(t),  may 
be  related  by  means  of  the  convolution  integral. 


2<t)  = 


H(t)  x(t  -  T)  dT 


(A.  3) 


where  T  is  the  variable  of  integration.   Squaring  Eq .  A. 3  yields 


x2(t)  = 

L  -oo  _i   <-  -oo 

Substituting  Eq.  A. 4  into  Eq .  A. 2  gives, 


H(T1)    z(t    -   T1)dT1 


e2(t)    =  x2(t)    -   2 


H(T)     z(t    -    T)dT    X(t) 


H(t)    z(t    -    T)dT 


H(T1)     Z(t     -     T1)dT1 


(A.  4) 


(A. 5) 
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By  definition,  the  mean  square  error  is, 

T 


2^   A    i-   -i- 
e  (t)   =    lim  2T 


e  (t)dt 


(A. 6) 


-T 


Substituting  Eq.   A. 5    into  Eq .    A. 6  yields: 


? 
e    (t) 


V£      2T 

X    oo 


x    (t)dt  - 


2  lim     - 

T->oo         z 


dt  H(T)z(t    -    T)x(t)dT 


+  lim  -  dt 

T-*co     2T 


H(T)z(t    -    T)dt  H(T1)z(t     -    T1)dT1 


(A. 7) 


The  definition  of   the   input  autocorrelation   function   is 

T 


cp      (t)    =   lim 
zz  T-»t» 


i     f 

2T     J 
-T 


z(t)z(t  +  T)dt 


(A. 8) 


Similarly,  the  crosscorrelation  function  between  the  input  and  the 
output  is 


'^(T>  '   itS?  21  [   z(t)5(t  +  T)dt 


(A.  9) 


and  the  autocorrelation  function  of  the  output  is, 


^(T)  -  li4m   2T  /  *(t)*(t  +  T)dt 

-T 


(A. 10) 


Now,  making  use  of  the  correlation  functions,  and  interchanging  the 
order  of  integration  in  Eq.  A. 7  yields, 

00 

e2(t)    =  cp      CO)    -   2     /     H(t)  cp      (T)dT  + 

Txx  Tzx 


00  00 

J  H(T)dT         I    H(Tl)    Cpzz(T     -    Tl)dT1 


(A. 11) 
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To  determine  the  system  function  that  minimizes  the  mean  square 
error,  assume  that  a  solution  exists  and  denote  this  solution  as 
H  (t).   Now  construct  a  weighting  function: 


m 


H(t)  =  H  (t)  +  6H.(t)  (A. 12) 


m 


l6 


where  H  (t)  is  the  assumed  solution,  H.(t)  is  any  arbitrary  realizable 

weighting  function,  and  6  is  a  parameter  which  may  be  varied  to  test 

whether  H  (t)  is  the  solution.   The  construction  is  such  that  at  6 
m 

equal  to  zero,  the  derivative  of  mean  square  error  with  respect  to  6 

must  be  equal  to  zero.   By  setting  this  derivative  equal  to  zero  for 

6  equal  to  zero,  a  condition  for  H  (t)  which  must  be  satisfied  in 

m 

order   for   it    to  be  a   solutionis   specified. 

Substituting  H(t)    as   given  by  Eq .   A . 12   into  Eq.   A. 11  and  differenti- 
ating with  respect   to  6   yields: 

00 

00  00 


0^J    =-    2    [h 


6(T)V(T)dT    +  Hm(T)dT  H    (T    >p  JfT    -    T    )* 


d6  J      6        Tzx  I       m  /       6      1  Yzz  1        1 


+ 

-00 


J     H.(T)dT  I      H     (T     )     Cp        (T     -    T     )dT 

Jo  Jmlzz  11 


03 


+  26    /    H6(T)dT     I     H6(t)  cpzz(T    -  T1)dT1  (A. 13) 

-00  -00 


Because  of  the  even  property  of  autocorrelation  functions  of  stationary 
signals , 

V(T1  "  T)  =<fzz(T-  V  (A-l4) 

which  leads  to 
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/". 


(T)dT  HgCrp    cpzz(T    -    T]_)dT1 


H6(T)dT     W^  -  Tl)dTl 


(A.  15) 


Substituting  Eq .  A. 15  in  Eq.  A . 13  and  setting 


e  (t)J  . 


d6 


=  0    at  6  =  0 


yields : 


2    H6(T)dT 


H  (T  )cp   (T  -  T  )dT   -  Cp   (T) 

m   1   zz      1    1    zx 


l.      -CO 


=  0     (A.  16) 


Since  H.  (T)  is  a  realizable  weighting  function,  it  must  be  zero  for 
values  of  t  less  than  zero.   The  only  way  that  Eq .  A. 16  can  be 
satisfied  for  T  ^  0  is  that  the  factor  in  the  brackets  be  equal  to 
zero.   Thus, 


H   (T,)'CD   (T  -  TjdT,  =  Cp   (T),     for  T  =  0 

m   1   zz      1    1   Tzx 


(A.  17) 


Equation  A. 17  is  the  famous  Wiener-Hopf  integral  equation. 

Strictly  speaking,  the  H  (t)  that  satisfies  Eq.  A. 17  produces  a 
stationary  value  of  the  mean  square  error  and  does  not  necessarily 
ensure  a  minimum.   It  can  be  shown  that  the  solution  of  Eq.  A. 17  does 
yield  a  minimum,  by  considering  the  second  derivative  of  the  mean 
square  error  with  respect  to  6.   Differentiating  both  sides  of  Eq.  A. 13 
yields 


lilt 


d6 


-2   /   H6(T)dT       H  (T  )  CPzz(T  -  T  )dT 
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This  is  the  same  as  twice  the  mean  square  value  of  the  noisy  input 
signal  filtered  by  a  weighting  function,  EL  (t).   Since  this  mean 
square  signal  can  never  be  negative,  this  shows  that  the  second 
derivative  of  the  mean  square  error  is  always  positive,  hence  the 
solution  of  Eq.  A . 17  corresponds  to  a  minimum. 
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APPENDIX  B 
DEVELOPMENT  OF  THE  BODE- SHANNON  SOLUTION 

Bode  and  Shannon  used  frequency  domain  techniques  to  solve  the 
filter  problem  by  converting  the  actual  signal  power  spectral  density 
to  that  of  white  noise,  and  then  operating  on  this  white  noise  signal. 
The  solution  may  be  described  by  considering  Figure  1;  x(t)   is  the 
signal,  v(t)   is  white  noise,  x(t)   is  the  estimate  of  the  signal,  and 
e(t)    is  the  error.   The  assumptions  are  that,  (1)  input  signal  and 
noise  are  uncorrelated  and  stationary,  and  (2)  all  time  functions  are 
Fourier  transformable. 

The  impulse  response  of  the  filter,  h(t),  is  chosen  on  the  basis 
of  minimum  mean  square  error.   The  error  signal  is  defined  as 


e(t)  =  x(t)  -  x(t) 
and  its  Fourier  transform  is, 


(B.l) 


E(UJ)  =  X(ou)  -  X(uo) 

=  fx(uu)  +  V(w)j   h(uo)  -  X(U)) 
=  fh(u))  -  lj  X(uo)  +  h(0J)V(uj) 


(B.2) 


Since  the  signal  and  noise  are  uncorrelated,  the  power  spectral  density 
of  the  error  is 


W  (u>)  =  ~  |h(uo)  -  1  I  2  W  (u>)  +  Ih(OJ)  I  2  W  (u>) 

e      2rr    '         '    x  '    v 


(B.3) 


Hence  the  mean  square  error  is 


e2(t)   -  i; 


h(0J)  -  l|2  W  (P)   +   |h(u))|2  W  (ui)    duJ   (B.4) 


This  is  to  be  minimized  by  the  proper  choice  of  h(uo) 
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Note  that  to  obtain  a  solution,  the  power  spectral  densities 
W  (UJ)  and  W  (UJ)  must  be  specified.   Let 


v 


h(o>)  =  C(uJ)e 


j6(uj) 


=  C(uj)   (CosQ(uj)  +  j  sin9(uJ)) 


(B.5) 


Substituting  into  Eq.  B.l  and  simplifying  yields, 


^-h 


C  (uu)W  (u>)  + 


C  (u>)  +  1  -  2C(uj)cos(0(uJ)) 


W  (uu)dUJ 


(B.6) 


Since   C(uj),   W    (uj) ,    and  W    (o>)    are  nonnegative   this   expression   is 
x  v 

minimized  when  cos   9(uo)   has    its    largest  value;    that   is,   when  0(uJ) 
equals    zero.      Thus, 


e2oPt(t)  =  h 


C    (tt)) 


W    (»)    +  W    (uj)        -    2C(u))W    (tu)   +W    (au) 

V  X  I  XX 


duo 
B.7) 


In  order   to   find   the  C(oj)    that  minimizes    this   expression,    complete 
the   square   of  the    terms    in  the  bracket.      Thus, 


e   opt(t)    : 


2tt 


C(ou)  VH    (ou)   +  W    (au)    - 

V  X 


W    (uj) 

X 


Vw    (uj)    +  W    (UJ) 

V  X 


W    (UJ)W    (tt)) 

XV 

W    (UJ)    +  W    (UJ) 

X  V 


duo 


Eq.    B.8   is  minimized  when 
C(UJ)    = 


W    (uj) 
x 


W    (uj)    +  W    (uj) 

X  V 


(B.8) 


(B.9) 
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with  minimum  mean  square  error  given  by 


CO 


W  (uu)W  (ou) 


e2  .  (t)  =  -J-    —^ 2 dU)  (Ba0) 

X         V 


Since  0(co)  equals  zero, 

W  (to) 

h(U))  =  C(CU)  =     - (B.ll) 

W  (u>)  +  W  (oo) 

X         V 

The  impulse  response  is  therefore  given  by 

CO 

f     W  (u>) 

h<t)  =  7    I    wx(°")X+wv(.)  ^  *         (B-L2) 

-co 

which,  in  general,  does  not  vanish  for  t  <  0.   Thus  Eq.  B.ll  is 
generally  not  physically  realizable. 

The  minimum  phase  filter  which  converts  a  white  noise  input  to  a 
spectral  density  equal  to  that  of  the  signal  plus  noise,  must  have 
a  transfer  function  of  magnitude 


B(U))  =  Vw  (oj)  +  W  (u)) 

X  V 

Since  Eq .  B.ll  gives  the  best  noncausal  operation  on  the  signal  plus 
noise,  the  best  noncausal  operation  on  white  noise  inputs  requires  that, 


W  (co)  \J   W  (u>)  +  W  (u)) 

V°>  "  -" wlS)  +  *%.)  <B-13> 

X  V 

This  however,  is  still  noncausal.   The  causal  filter  is  then  constructed 
as  follows: 

0         ,  t  <  0 
gjt)   = 

h(ju(t)      ,  t  £  0  (B.14) 
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The  desired  causal  filter  for  operation  on  signal  plus  noise 
inputs  then  has  a  transfer  function  given  by: 

G  (UJ) 

B(U°      "       Vw  (u>)  +  W  (ou) 

X         V 
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APPENDIX  C 


POWER  SPECTRAL  DENSITY  PROGRAM 
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POWER  SPBCTRAL  DENSITY  PROGRAM 
FORTRAN  IV  G  LEVEL  0,  MOD  0  MAIN  DATE  =  67328  18/15/48 


0CC1 


c   ~p ml! Ablflem<libE/iF2 Input  array  to  harm  must 

C   BE  A  POWER  OF  TkO 


T(C2  CCMPLEX*8  A<2048,1,1» 

0003  DIMENSION  S< 1024  J , INV( 1024) 

0004  DIMENSION  M  (  3  I  ,  X  ( <*  1  00  I  ,  EU  200  )  ,  HX  (  200  )  ,  T  I  MEX  <  200  I 

0005  KK=11 
0CC6  T=C01 

0CC7  PIT=6.2b31853 

00C8  QD=1.0/T 

0Cu9  DEV=SQRT<QC) 

0010  R=1.0/T 

0011  OEX=SCRT(R) 
l»012  PM  =  EXP(-T) 

0013  GAPMA=1.0-PHI 

0014  IX=4C9734751 

0015  CALL  GAUSSIIX,DEV,0.0,W) 

0016  x(l)=GAMMA*h 

0017  M<1)=KK 
CO  18  M(2)=0 
Col9  M<3)=0 

)C2C  M  =  2**KK 

0021  N2=l 

00  2  2  N3=l 

0C2  3  DC  2  U=1,M 

0(24  DC  2  I2=UN2 

0025  DC    ?    13=1, N3 

OC 2  6  2  A<  I  1,I2,I3)  =  (0.C ,0.0) 
C    SIGNAL  IS  NOVv  GENERATED 

0027  DC  3  11=1, M 

0028  00  3  12=1, K? 
'029  Of  3  13  =  1,  N3 

003C  CALL  GAUSS(  I  X  ,0E V , 0  .0 , W ) 

0031  XI  ll*l)=PHI*X( I1)*W*GAMMA 

UC32  3  A(  I1,12,I3)  =  X( II) 

f.   CALCULATE  MFAN  AND  MEAN  SQUARE  OF  SIGNAL 

00  3  3  AX=0.0 

0034  BZ  =  C0 

0C35  DC  50  1  =  1,  M 

0036  5C  AX=AX+X(  I  ) 

0037  XBAR  =  AX/M 

0038  DO  51  1=1, Nl 

0039  51  BZ=BZ+( X< I )-XBAR)**2 

0040  XKS  =  BZ/M 

0041  WRITF(6,92)XBAR 

0042  92  F0RMATIT6,  'SIGNAL  AVERAGE  VALUE  =  '♦T32,E13.6) 

0043  WRITE(6,95)XMS 

0~44  95  F0RMAKT6, 'MFAN  SQUARE  SIGNAL  =  «,T30,E13.6) 

C   OBTAIN  THE  FAST  FOURIER  TRANSFORM 

C045  CALl    HAftM*A,M,INV,S,l,IP6R«) 

CC6  NX=35 

C   NORMALIZE  AND  PREPARE  INPUTS  FOR  PLOTTING 

0047  DC  42  1=1, NX 

0048  B( I )=CA8S( A( 1,1,1) )**2 
0G49  BX(I )  =  (B( I  )*2.0)/(Nl*QD) 
0O5C  TIMEXI  I)  =  ( (  I-1)*PIT)/<N1*T) 

0051  HRITE(6,63)TIMEX{ I  ),BX( I ) 

0052  63  FCRMAl(T6rE13^6,T2C,E13.6) 

0053  42  CONTINUE 

0054  READ(5,10  5)(ITITLE!I),I=1,12) 

0055  105  FCRMAT(6A8) 

0056  CALL  DRAW(NX,TIMEX,BX,0,0» •     • , I T I TLE ,0.0,0.  OtO,OtO^O, 8,0» 1 ,L A ) 
CC57  END 
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APPENDIX  D 


WIENER-KOLMOGOROV  PROGRAM 
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WIENER-XOLMOGOROV    FILTER   PROGRAM 

FCklRAN  IV    C    LEVEL    Oi    MGD    C                                 MAIN                                      DATE    =   67330                         14/00/57 

88c  2  If  4fe£§  I  hi l  JV^  lo^lIXMAT(3100>fi«3l004rZfMO»»TTIMEI31001 

CCC3  OIMENSICN    EMS ( 3  ICO) , VXX( 310  I ,ES< 3100  I 

COCA  N=3CCC 

0CC5  NA=N+5 

COCfc  IX=ACS73A751                     —            — 

0<X7  T  =  0.C1 

COCP  QC=1.C/T 

COi<;  OEV=SCRT<GC) 

001C  RM.O/T 

CC11  DEX=SCRT(R) 

0012  EMS(l)=C.O 

0C13  PHI=FXP<-T) 

G01A  GAHMA=l.O-PHl 

CC15  RCT=SCRT(2.C) 

0C16  XPC=RCT*T 

0017  PHK=EXP(-XFC) 

0018  GAMg=l„Q-PRK  -        _— _.                        

C01c;  V»RITE(6,2C) 

00 2C  20    FORMAT < T6  ,  •  K  •  ,  T16 , • X ( K ) • ,T36 , 'XHATI K ) • , T 58 , • Z( K » • ,T80, • EMS(K ) • ) 

0021  XFAT(1>=0.0 

0022  CALL  GAUSS(IX,DEV,C.0,WI 

0023  XI  1)=GAMMA*V< 

C      GENERATE  SIGNAL 

002A  DC  10  K=1,NA 

0025     CALL  GAUSS!  I)1  ,OEV,0»0,M) — 

CG26  X(K+1)=PHI*X(K)*V»*GAMMA 

C027  10  CONTINUE 

C      ACD  MEASUREMENT  NOISE  AND 

C      CCMPUTE  ESTIMATE  Of  SIGNAL 

002P  DC  11  K=1,NA 

CC29  CALL  GAUSS! I  X ,CE X ,0 .0, V) 

CC3C  Z(K)=X(K)*V 

0031 XRAT+K*l)=P*K*XfcO-T4*UGAMB«ZU)*(  U0^42.O*R0T»  ) 

C032  E<K)=XHAT(K)-X(K) 

CC33  11  CCNTINUE 

CC34  KT=2 

C      CCMPUTE  MEAN  AND  MEAN  SQUARE  ERROR 

CO  3 5  DO  8  5  K=1,N 

L'036  YA  =  C.O 

OC27  Y6=C.O 

0038  DC  25  1=1, KT        

CC?<;  YA=YA+F(I> 

CGAC  25  CCNTINUE 

C?A1  YPAR=YA/KT 

CCA2  DC  26  1=1, KT               

0CA3  Ye=YP+( E(  I  l-YBAR 1**2 

CC<<A  26  CONTINUE 

COA1-  EMS(K+1)  =  YB/KT 

00 A 6  K4-KT+1 


CCA7  V<PITE(6f30)K,X(K),XHAT(K)1Z(K)  ,EMS(K) 

CJAE  30  FCRMAT(2X,!A,5X,E13.6,9X,E13.6,9X,E13.6,9X,F13.6) 

OCA'S  e5  CCNTINUF 

C  PREPARE  ITEMS  FOR  PLOTTING 

CC5C  CC  7C  1  =  1, N 

C051  TIME(I)=I-1 

CC52  7C  CCNTINUE 

C"53  -OC-HO  A=irr*i-                   

OOtA  110  ESII  )  =  E(  I)**2 

C.055  KC=1 

CL56  NPT=3C0 

f057  DO  A3  I=l,NA,10 

CC58  YXX(KQ)=EMS( I ) 

C'5<;  A3  KC  =  KQ+1 

<-.,6C  REACI5.1C1  )(  ITITLEI  I  J, 1  =  1, 12) 

0061  1G1  FCRVATl6Ag) 


)C62  CAIL  rRAW(NPT,TIME,YXX,0,C,«     • , I TI TLE ,0.0 ,0 .0,Of 0 ,0,0 , 8, 9 , 1 , L A ) 

0Lf3  ENO 
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APPENDIX  E 


KALMAN  PROGRAM 
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KALMAN  FILTER  PROGRAM 


FORTRAN 

m\ 

00C3 

00C4 

00C5 

0CC6 

0CC7.._ 

0CC8 

0CC9 

0010 

0C11 

0012 

0n  13 

0014 

0015 

0016 

00  17 

p.ni  a 

CC19 

CC2C. 

Or  ?1 

0022 

C>2^ 

0024 

oC  2  5 

0026 

0C27 

0028 

0C29 

0C3C 

00?  1 

0132 

00  3  3 

0C34 

0  35 

0.  36 

O037 

oo  3  n 

0C39 

0:40 

C  '•  4 1 

c ■ :  ^  2 

0C43 

0044 

0045 

0C46 

(..  C  4  7 

V  4  8 

0049 

00  5  C 

0C51 

0i.52 

0053 

0  )54 

-j  •  (,  c. 

j   '.-  u 

0057 

v.  .CP 

00  59 

;.  0  6  C 

006  1 

0062 

006  3 

utt4 

^  <  c. 

a  'i '  - 

IV-6   L§V*4-G-r    MOD    0 MA4AI 0AT6    -   6  7*31— — -         -     14/20/4* 


mm  kHMHU 


DIMENSION  xniOO  ).  XHAT<  31001  «E<  3  100  >,Z(  3100)  ,T  I  ME  I  3100  1  ,  EMS<  3100  I 
DIMENSION  YXX(310)TGX(3IOT»PM3100),P(3ir>0),G<  31001, ES(3100I 

T=0.01  ~ 

N=3000 
1*A=N*5    .     .  ..   

Q0=1.0/T 

DEV*SCRT(QC) 

R-l.O/T ♦- 

DEX=SGRT(R) 

IX=4C9734751  

EHSM)  =  P.O  

PHI=EXP(-T) 
_    GAKMA=1.0-PHI  

XHAT(  l)=0.C 

Q=GAMMA*GAM>'A*<DEV**2) 

PM*  =  10.0 

NP=NA+2 
C  GENERATE    GAIN    SEQUENCE 

DV~*i~?    K=1,NR 

G(K)=P(K)/(P(K )«-R ) 

PKIK)=(1.C-G(K))*P(K1         

57    P(K*1)=PHI*PM*PK(K)+Q 

WPITEI6.2C ) 
20    F€RMAT«T6,«K«  ,  T  1  3  ,  •  X4  K  *  •  TT32~r«  XHAT4K  )  «,T50,«7<  K  )  •  ,  T70rl€MS  <K  »  •  ,  T90 
1,«G(K)« I 

CALL    GAUSS< I  X  ,DE  V  ,  0.0  ,  I*  ) 

XtTT=GAMMA*V 
C      GENERATE  SIGNAL 

CC  1C  K=1,NA 

CALL  GAUSSU  X,OEV,0.C,W) 

X(K+  1  )  =  PHI*X (K  )  +  U*GAMMA 
-   10  CCNTINUF 

CALL    GAUSS( I  V, R6X, 0.0, V) 

Z(  n  =  x<  i)  +  v 

C  ATD    MEASUPE^FNT    NOISE 

C      CCMPUTF  ESTINATF.  OF  SIGNAL 

DC  11  K=1,NA  

CALL  GAUSS!  1 X,DFX,i .O.V) 

Z(K+1  (  =  X(K4l  )+V 

XHAT<K«-H  =  PHI*XHAT<K)+G»K*l>*<^K«-l  >-PHI*XHAT(  K)  ) 

E(K)=XHAT(K)-X(K) 
11  CONTINUE 

KT=2 
C      CALCULATE:  f*EAM  AND  MEAN  SQUARE  OF  ERROR 

CC  85  K=) ,N 

YA^.C 

YR=  ;.c 

CC  25  1=1, KT 

YA  =  YA  +  E(  I  ) 
2  5  CONTINUE 

YPAR=YA/KT 

DO  26  1  =  1, KT 

YB=YH+(£{ I  )-YBAR)**2 
26  CTNTIN'UF 

EVS(K+1 )=YP/KT 

KT=KT+1 

WRITE(6,_i~)K,X<K),XHAT(K),Z(K),EMS(K),G(K) 
30  FCRMATCX,I4,5X,E13.6t9X,E13.6,9X,E13.6,9X,F13.6,9X,E13.6) 
F5  CCNTINUE 
C     PRFPARE  ITEMS  FCK  PLOTTING 

DC  7C  1=1, N 

TIME!  I  1  =  1—1 
70  CCNTINUE 

DO  110  1  =  1  »N 
110  ESI  I  )  =  E( I  )**2 

KQ=1 

NPT=3GC 

DO  43  I=1,NA,1C 

YXX(KQ)=Ef S( I ) 
43  KC"  =  KC+1 

READ (5, 10  1  )(  IT  I  TIF (  I  )  ,1  =  1  ,12) 
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KALMAN  FILTER  PROGRAM  

FCKTU&K  IU   C-  LEV6L -0-*-MGO    0     -        -•-      -WW DAT€«67*34 1-6/80/48 

%Ybl.  I01  !:XfL*^Sft?ilpT,TiHefyxxfOie-f*         ,»!TtTLe,o.o,o.o,o,0te,o,8,9,itLA> 

^J7c  cn~4«    I=1,NA,1C 

C071  GX(KL»=G(!)  

C  )72  <«<»    KL  =  KL+1 

0)74  CALL    DRAhfKPT,TlWEtGx7CrfOT»            *7KT!TLE  ,0.0,0.  0,0,0,0, 0, 8 ,9, 1  ,L  A ) 

C075  END 
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APPENDIX  F 


NO  FILTER  PROGRAM 
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NO    FILTER    PROGRAM 

FCPTKAN     IV    G    LEVEL    Ct    MOO    0  MAIN                                       DATE    =   67330                          13/54/55 

88£  \  ol^E-NilON  l  \\\  I  H\    XHATf  51f0  \  ,  E-m^-fT*<  3100  >  ,T  fM€t  31001 

0013  DIMENSION    EKS( 31C0 ) , YXXI 310) , WI 3100  J , VI 3 100 ) ,E S ( 3100 ) 

0CC4  N=3000 

0CC5  NA=N*5 

00T6  T=0.01 

0017  QC=i.O/T 

OCfe  OEV=SQRTICC) 

C0C9  R=1.0/T 

C01C  — —  C€X*SQ«mU                              — ^— 

0011  IX=4C9734751 

oci2  xm=o.g 

0C12  EMSI1)=G.0 

0014  PHI  =  EXPI-T» 

0015  GAMMA=l.C-PHl 
GC16  WRITEI6.2C) 

0017  20  FORMAT! T6  ,  •«  •  ,  T 14  ,  • XI  K » • , T 34 , • 2  J K» • tT54, • EMS  IK ) • , T74  ,  «WI K  )  • , T94, 

1«V4K)M  

GC1P  CC  11  K=1,NA 

0019  11  CALL  GAUSSI IX.DEV,C.O,WIK)) 

C  GENERATE  SICNAL 

C  ACD  MEASUREMENT  NOISE  AND 

C  COMPUTE  ESTIMATF  OF  SIGNAL 

0C2C  CO  10  K=1,N 

0021  XIK*1)  =  PHI*XIK)*MK*1»*GAMMA 

0022  ZIK  )  =  VIK1*XIK j 
00 2 A  EIK)=Z(K>-X{K) 

0025  1C  CONTINUE 

C  COMPUTE  MEAN  AND  MEAN  SQUARE  ERROR 

0026  KT=2 

OC27  CC  85  K=l ,N 

0G2P  YA=O.C 

0C29 ¥«=-Q,£ — — 

0U3C  DC  25  1=1, KT 

0C21  YA=YA»EII> 

0032  25  CONTINUE 

0CJ3  YBA«=YA/KT 

C034  CC  26  1=1, KT 

0035  YB=VR*( EI  I  )-YBAR)**2 

0036  26  CONTINUE 

0037  €MS<K*1)=YB/KT 
CC3B  KT=KT*1 

CC39  taRITE(6,3ClK,X(K),Z<K),EMS(K)tW(KI,V(K> 

004  -  30  FrRMAT(2X,l4,4X,E13.6,7X,E13.6,7X,E13.6,7X,E13.6,7X,E13.6J 

0041  85  CONTINUE 

C  PREPARE  HEMf  FOR  PLOTTING 

ol42  CC  7C  1  =  1,  N 

CCA?  TIMtm  =  I-1 

0044  7G  CCMINUE 

C>04<:-  CC  110  1  =  1,  N 

004(  110  ESI  I  )  =  E( I »**2 

0047  KC=1 

0048  NPT=300 

0049  DC  43  I  =  1,NA,10 
005G  YXXIKQ»=EfS< I) 
0051  43  KC=KC+1 

0"b2  REAOI 5,1011 1  IT  1TLE  U  >  ,  I  =  1 ,  12  » 

0(>53  101    FCPMAU6A8) 

0"54  CALL    DRAW(NPT,TIME,YXX,0,C,«             •  ,1 T  I  TLE  ,0.0  ,0 .0 ,0 ,0 , 0,0 , 8, 9, 1 , L A > 

Ci'1"5  END 
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